// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtStatisticalTrendsFilter.C
// ************************************************************************* //

#include <avtStatisticalTrendsFilter.h>
#include <ImproperUseException.h>
#include <vtkDataSet.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
#include <vtkCellType.h>
#include <vtkIdList.h>
#include <vtkPoints.h>
#include <vtkDataArray.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>

// Use of named selection
#include <avtCallback.h>
#include <avtContract.h>
#include <avtParallel.h>
#include <avtIdentifierSelection.h>

#include <map>
#include <string>

// ****************************************************************************
//  Method: avtStatisticalTrendsFilter constructor
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Fri Jan 25 11:02:55 PDT 2008
//
// ****************************************************************************

avtStatisticalTrendsFilter::avtStatisticalTrendsFilter() :
  validTimeAxis(true), lastTimeAxisValue(-1.0e12)
{
    haveData = true;
    sumX_ds      = 0;
    sumY_ds      = 0;
    sumX2_ds     = 0;
    sumY2_ds     = 0;
    sumXY_ds     = 0;
    slope_ds     = 0;
    intercept_ds = 0;
    out_ds       = 0;

    numTypes = 6;

    strcpy( typeString[0], "Sum" );
    strcpy( typeString[1], "Mean" );
    strcpy( typeString[2], "Variance" );
    strcpy( typeString[3], "Standard Deviation" );
    strcpy( typeString[4], "Slope" );
    strcpy( typeString[5], "Residuals" );
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter destructor
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Fri Jan 25 11:02:55 PDT 2008
//
//  Modifications:
//
// ****************************************************************************

avtStatisticalTrendsFilter::~avtStatisticalTrendsFilter()
{
  if( trend_ds.size() )
  {
    for(size_t i=0; i<trend_ds.size(); ++i )
      trend_ds[i]->Delete();

    trend_ds.clear();
  }
}


// ****************************************************************************
//  Method:  avtStatisticalTrendsFilter::Create
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Fri Jan 25 11:02:55 PDT 2008
//
// ****************************************************************************

avtFilter *
avtStatisticalTrendsFilter::Create()
{
    return new avtStatisticalTrendsFilter();
}


// ****************************************************************************
//  Method:      avtStatisticalTrendsFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Fri Jan 25 11:02:55 PDT 2008
//
// ****************************************************************************

void
avtStatisticalTrendsFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const StatisticalTrendsAttributes*) a;
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::Equivalent
//
//  Purpose: Returns true if creating a new
//      avtStatisticalTrendsFilter with the given parameters would
//      result in an equivalent avtStatisticalTrendsFilter.
//
//  Programmer: childs -- generated by xml2avt
//  Creation:   Fri Jan 25 11:02:55 PDT 2008
//
// ****************************************************************************

bool
avtStatisticalTrendsFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(StatisticalTrendsAttributes*)a);
}

// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::ExamineContact
//
//  Purpose: Examine the contract to get the current state of the time
//    slider. The time slider state is needed in case that the start
//    and end time are relative to the time slider.
//
//  Programmer: Oliver Ruebel
//  Creation:   May 07, 2009
//
//    Oliver Ruebel, Thu May 11 10:50
//
// ****************************************************************************
void
avtStatisticalTrendsFilter::ExamineContract(avtContract_p in_contract)
{
    // Call the examine contract function of the super classes first
    avtPluginFilter::ExamineContract(in_contract);
    avtTimeLoopFilter::ExamineContract(in_contract);
    avtDatasetToDatasetFilter::ExamineContract(in_contract);
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::ModifyContract
//
//  Purpose:
//      Creates a contract the removes the operator-created-expression.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
// ****************************************************************************

avtContract_p
avtStatisticalTrendsFilter::ModifyContract(avtContract_p in_contract)
{
    avtDataRequest_p in_dr = in_contract->GetDataRequest();
    outVarName = in_dr->GetOriginalVariable();

    in_dr->SetUsesAllDomains(true);

    if( strncmp(outVarName.c_str(), "operators/StatisticalTrends/",
                strlen("operators/StatisticalTrends/")) == 0)
    {
      std::string operatorWithVar =
        outVarName.substr(strlen("operators/StatisticalTrends/"));

      for (int t = 0; t < numTypes; ++t)
      { 
        if( strncmp(operatorWithVar.c_str(), typeString[t],
                    strlen(typeString[t])) == 0)
        {
          if( t != atts.GetStatisticType() )
            atts.SetStatisticType( (StatisticalTrendsAttributes::StatisticTypeEnum) t );

          std::string justTheVar =
            operatorWithVar.substr(strlen(typeString[t])+1);

          outVarName = justTheVar;

          avtDataRequest_p out_dr =
            new avtDataRequest(in_dr, justTheVar.c_str());

          return avtPluginFilter::ModifyContract( new avtContract(in_contract, out_dr) );
        }
      }
    }

    return avtPluginFilter::ModifyContract(in_contract);
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Tells output that we have a new variable.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
// ****************************************************************************

void
avtStatisticalTrendsFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes &in_atts = GetInput()->GetInfo().GetAttributes();
    avtDataAttributes &out_atts = GetOutput()->GetInfo().GetAttributes();

    if( outVarName != "" )
    {
      std::string fullVarName = outVarName + " " +
        std::string( typeString[(int) atts.GetStatisticType()] );

      out_atts.RemoveVariable(in_atts.GetVariableName());
      
      if( !out_atts.ValidVariable(fullVarName) )
      {
        out_atts.AddVariable((fullVarName).c_str());
        out_atts.SetActiveVariable(fullVarName.c_str());
        out_atts.SetVariableDimension(1);
        
        out_atts.SetVariableType(AVT_SCALAR_VAR);
      }
    }

    avtPluginFilter::UpdateDataObjectInfo();
    avtTimeLoopFilter::UpdateDataObjectInfo();
    avtDatasetToDatasetFilter::UpdateDataObjectInfo();
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::InitializeTimeLoop
//
//  Purpose: Set the start, stop, and strides.
//
//  Programmer: Allen R. Sanderson
//  Creation:   May 07, 2011
//
// ****************************************************************************
void
avtStatisticalTrendsFilter::InitializeTimeLoop(void)
{
    if( atts.GetStatisticType() == StatisticalTrendsAttributes::Sum)
      SetNumberOfIterations(1);
    else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Mean)
      SetNumberOfIterations(1);
    else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Variance)
      SetNumberOfIterations(2);
    else if( atts.GetStatisticType() == StatisticalTrendsAttributes::StandardDeviation)
      SetNumberOfIterations(2);
    else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Slope)
      SetNumberOfIterations(1);
    else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Residuals)
      SetNumberOfIterations(2);

    int numStates = GetInput()->GetInfo().GetAttributes().GetNumStates(); 
    int startTimeSlice;
    int stopTimeSlice;

    if( atts.GetStartIndex() < 0 )
    {
      std::string msg(GetType());
      msg = msg + ": Start index/Number of slices must be positive.";
      EXCEPTION1(ImproperUseException, msg);
    }

    if( atts.GetStopIndex() < 0 )
    {
      std::string msg(GetType());
      msg = msg + ": Stop index/Number of slices must be positive.";
      EXCEPTION1(ImproperUseException, msg);
    }

    // If the trend type is absolute.
    if( atts.GetStartTrendType() == StatisticalTrendsAttributes::Absolute )
    {
        startTimeSlice = atts.GetStartIndex();
    }

    // If the trend type is relative to the current time slider state
    else //if( atts.GetStartTrendType() == StatisticalTrendsAttributes::Relative )
    {
        int startOffset = atts.GetStartIndex();
        startTimeSlice = currentTime - startOffset;
    }

    // If the trend type is absolute.
    if( atts.GetStopTrendType() == StatisticalTrendsAttributes::Absolute )
    {
         stopTimeSlice = atts.GetStopIndex();
    }

    // If the trend type is relative to the current time slider state
    else //if( atts.GetStopTrendType()  == StatisticalTrendsAttributes::Relative )
    {
        int stopOffset = atts.GetStopIndex();
        stopTimeSlice = currentTime + stopOffset;
    }

    // Check if the times are valid and correct if needed
    if (startTimeSlice < 0)
      startTimeSlice = 0;
    
    if( startTimeSlice >= numStates )
      startTimeSlice = numStates-1;
    
    if( stopTimeSlice < 0 )
      stopTimeSlice = 0;
    
    if (stopTimeSlice >= numStates)
      stopTimeSlice = numStates-1;

    // Update the start and end frame as well as the stride.
    SetStartFrame(startTimeSlice);
    SetEndFrame(stopTimeSlice);
    SetStride( atts.GetStride() );

    // Clean-up the current output dataset if necessary
    if( trend_ds.size() )
    {
      for(size_t i=0; i<trend_ds.size(); ++i )
        trend_ds[i]->Delete();
      
      trend_ds.clear();
    }

    avtTimeLoopFilter::InitializeTimeLoop();
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::Execute
//
//  Purpose: Defines what it means for this filter to "Execute". This
//    is where the actual iteration over time happens. This functions
//    is overwritten here to allow the dynamic setting of start and
//    end-time of the iteration.  The iteration over time is performed
//    using avtExecuteThenTimeLoopFilter::Execute(void)
//
//  Programmer: Allen R. Sanderson
//  Creation:   May 07, 2011
//
// ****************************************************************************

void
avtStatisticalTrendsFilter::Execute(void)
{
    avtDataTree_p tree = GetInputDataTree();

    // Ask for the dataset
    int nds;
    vtkDataSet **dsets = tree->GetAllLeaves(nds);

    if (nds == 0)
    {
        // Free the memory from the GetAllLeaves function call.
        delete [] dsets;

        return;
    }

    int nds2 = nds;
    SumIntAcrossAllProcessors(nds2);

    if (nds2 < 1 || 1 < nds)
    {
        // Free the memory from the GetAllLeaves function call.
        delete [] dsets;

        EXCEPTION1(ImproperUseException, "Filter expected only one vtkDataSet"
                                         " in avtDataTree");
    }

    vtkDataSet *curr_ds = dsets[0];

    // Free the memory from the GetAllLeaves function call.
    delete [] dsets;

    // Time axis value - only used when time sampling.
    double timeAxisValue;

    if( atts.GetTrendAxis() == StatisticalTrendsAttributes::Time &&
        GetInput()->GetInfo().GetAttributes().TimeIsAccurate() )
      timeAxisValue = GetInput()->GetInfo().GetAttributes().GetTime();
    
    else if( atts.GetTrendAxis() == StatisticalTrendsAttributes::Cycle &&
             GetInput()->GetInfo().GetAttributes().CycleIsAccurate() )
      timeAxisValue = GetInput()->GetInfo().GetAttributes().GetCycle();
    else
      timeAxisValue = currentTime;

    if( GetIteration() == 0 && lastTimeAxisValue >= timeAxisValue )
      validTimeAxis = false;

    // std::cerr << timeAxisValue << "  " << lastTimeAxisValue << "  "
    //        << validTimeAxis << "  " << currentTime << "  "
    //        << GetInput()->GetInfo().GetAttributes().GetTime() << "  "
    //        << GetInput()->GetInfo().GetAttributes().GetCycle() << "  "
    //        << std::endl;

    lastTimeAxisValue = timeAxisValue;

    // Sumation iteration
    if( GetIteration() == 0 )
    {
      // First time step of first iteration so make a copy. 
      if( currentTime == GetStartTime() )
      {
        if( atts.GetStatisticType() == StatisticalTrendsAttributes::Sum ||
            atts.GetStatisticType() == StatisticalTrendsAttributes::Mean )
        {
          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[0]->DeepCopy( curr_ds );
          sumY_ds = trend_ds[0];

          out_ds = trend_ds[0];
        }
        else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Variance ||
                 atts.GetStatisticType() == StatisticalTrendsAttributes::StandardDeviation )
        {
          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[0]->DeepCopy( curr_ds );       
          sumY_ds = trend_ds[0];

          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[1]->DeepCopy( curr_ds );
          out_ds = trend_ds[1];
        }

        else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Slope ||
                 atts.GetStatisticType() == StatisticalTrendsAttributes::Residuals )
        {
          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[0]->DeepCopy( curr_ds );
          sumX_ds = trend_ds[0];

          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[1]->DeepCopy( curr_ds );
          sumY_ds = trend_ds[1];

          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[2]->DeepCopy( curr_ds );
          sumX2_ds = trend_ds[2];

          trend_ds.push_back ( curr_ds->NewInstance() );
          trend_ds[3]->DeepCopy( curr_ds );
          sumXY_ds = trend_ds[3];

          if( atts.GetStatisticType() == StatisticalTrendsAttributes::Slope )
          {
            trend_ds.push_back ( curr_ds->NewInstance() );
            trend_ds[4]->DeepCopy( curr_ds );
            out_ds = trend_ds[4];
          }

          else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Residuals )
          {
            trend_ds.push_back ( curr_ds->NewInstance() );
            trend_ds[4]->DeepCopy( curr_ds );
            slope_ds = trend_ds[4];

            trend_ds.push_back ( curr_ds->NewInstance() );
            trend_ds[5]->DeepCopy( curr_ds );
            intercept_ds = trend_ds[5];

            trend_ds.push_back ( curr_ds->NewInstance() );
            trend_ds[6]->DeepCopy( curr_ds );
            out_ds = trend_ds[6];
          }
        }

        // Zero out the sums.
        if( curr_ds->GetPointData()->GetScalars() )
        {
          unsigned int tPoints = trend_ds[0]->GetNumberOfPoints();
          unsigned int nPoints = curr_ds->GetNumberOfPoints();
          
          if( nPoints && nPoints == tPoints )
          {
            double val = 0;
            
            // Traverse all point data
            for( unsigned int i=0; i<tPoints; ++i )
            {
              if( sumX_ds )
                sumX_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( sumY_ds )
                sumY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( sumX2_ds )
                sumX2_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( sumY2_ds )
                sumY2_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( sumXY_ds )
                sumXY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( slope_ds )
                slope_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( intercept_ds )
                intercept_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

              if( out_ds )
                out_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
            }
          }
        }

        if( curr_ds->GetCellData()->GetScalars() )
        {
          unsigned int tCells = trend_ds[0]->GetNumberOfCells();
          unsigned int nCells = curr_ds->GetNumberOfCells();
          
          if( nCells && nCells == tCells )
          {
            double val = 0;
            
            // Traverse all cell data
            for( unsigned int i=0; i<tCells; ++i )
            {
              if( sumX_ds )
                sumX_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( sumY_ds )
                sumY_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( sumX2_ds )
                sumX2_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( sumY2_ds )
                sumY2_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( sumXY_ds )
                sumXY_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( slope_ds )
                slope_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( intercept_ds )
                intercept_ds->GetCellData()->GetScalars()->SetTuple(i, &val);

              if( out_ds )
                out_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
            }
          }
        }
      }

      if( curr_ds->GetPointData()->GetScalars() )
      {
        int tPoints = trend_ds[0]->GetNumberOfPoints();
        int nPoints = curr_ds->GetNumberOfPoints();
        
        if( nPoints && nPoints == tPoints )
        {
          // Get the next point and update its coordinates if necessary
          vtkDataArray *curr_scalars = curr_ds->GetPointData()->GetScalars();

          double val;

          // Traverse all point data
          for( unsigned int i=0; i< (unsigned int)nPoints; ++i )
          {
            if( sumX_ds )
            {
              val = *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) +
                timeAxisValue;
              
              sumX_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumY_ds )
            {
              val = *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i)) +
                *(curr_scalars->GetTuple(i));
              
              sumY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumX2_ds )
            {
              val = *(sumX2_ds->GetPointData()->GetScalars()->GetTuple(i)) +
                timeAxisValue * timeAxisValue;
              
              sumX2_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumY2_ds )
            {
              val = *(sumY2_ds->GetPointData()->GetScalars()->GetTuple(i)) +
                *(curr_scalars->GetTuple(i)) * *(curr_scalars->GetTuple(i));
              
              sumY2_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumXY_ds )
            {
              val = *(sumXY_ds->GetPointData()->GetScalars()->GetTuple(i)) +
                timeAxisValue * *(curr_scalars->GetTuple(i));
              
              sumXY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
            }
          }
        }
      }
      
      if( curr_ds->GetCellData()->GetScalars() )
      {
        int tCells = trend_ds[0]->GetNumberOfCells();
        int nCells = curr_ds->GetNumberOfCells();
        
        if( nCells && nCells == tCells )
        {
          // Get the next point and update its coordinates if necessary
          vtkDataArray *curr_scalars = curr_ds->GetCellData()->GetScalars();

          double val;

          // Traverse all cell data
          for( unsigned int i=0; i< (unsigned int)nCells; ++i )
          {
            if( sumX_ds )
            {
              val = *(sumX_ds->GetCellData()->GetScalars()->GetTuple(i)) +
                timeAxisValue;
              
              sumX_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumY_ds )
            {
              val = *(sumY_ds->GetCellData()->GetScalars()->GetTuple(i)) +
                *(curr_scalars->GetTuple(i));
              
              sumY_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumX2_ds )
            {
              val = *(sumX2_ds->GetCellData()->GetScalars()->GetTuple(i)) +
                timeAxisValue * timeAxisValue;
              
              sumX2_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumY2_ds )
            {
              val = *(sumY2_ds->GetCellData()->GetScalars()->GetTuple(i)) +
                *(curr_scalars->GetTuple(i)) * *(curr_scalars->GetTuple(i));
              
              sumY2_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
            }
            if( sumXY_ds )
            {
              val = *(sumXY_ds->GetCellData()->GetScalars()->GetTuple(i)) +
                timeAxisValue * *(curr_scalars->GetTuple(i));
              
              sumXY_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
            }
          }
        }
      }
    }

    // Sum of squares iteration
    else if( GetIteration() == 1 )
    {
      // First time step of second iteration so get the averages.
      if( currentTime == GetStartTime() )
      {
        // Calculate the average.
        if( curr_ds->GetPointData()->GetScalars() )
        {
          unsigned int tPoints = trend_ds[0]->GetNumberOfPoints();
          unsigned int nPoints = curr_ds->GetNumberOfPoints();
          
          if( nPoints && nPoints == tPoints )
          {
            // Traverse all point data
            for( unsigned int i=0; i<tPoints; ++i )
            {         
              if( atts.GetStatisticType() ==
                  StatisticalTrendsAttributes::Variance ||
                  atts.GetStatisticType() ==
                  StatisticalTrendsAttributes::StandardDeviation )
              {
                double val =
                  *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i));
                val /= (double) GetNFrames();
              
                sumY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
              }
              else if( atts.GetStatisticType() ==
                  StatisticalTrendsAttributes::Residuals)
              {
                // Slope
                double val =
                  (*(sumXY_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                   (*(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                    *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i))) /
                   (double) (GetNFrames())) /
                  (*(sumX2_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                   (*(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                    *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i))) /
                   (double) (GetNFrames()));
                
                slope_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

                // Average values
                val = *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i)) /
                  (double) GetNFrames();
                sumY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

                val = *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) /
                  (double) GetNFrames();
                sumX_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

                // Intercept
                val= *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                  *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                  *(slope_ds->GetPointData()->GetScalars()->GetTuple(i));

                intercept_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
              }
            }
          }
        }

        if( curr_ds->GetCellData()->GetScalars() )
        {
          unsigned int tCells = trend_ds[0]->GetNumberOfCells();
          unsigned int nCells = curr_ds->GetNumberOfCells();
          
          if( nCells && nCells == tCells )
          {
            // Traverse all point data
            for( unsigned int i=0; i<tCells; ++i )
            {
              if( atts.GetStatisticType() ==
                  StatisticalTrendsAttributes::Variance ||
                  atts.GetStatisticType() ==
                  StatisticalTrendsAttributes::StandardDeviation )
              {
                double val =
                  *(sumY_ds->GetCellData()->GetScalars()->GetTuple(i));
              
                val /= (double) GetNFrames();
              
                sumY_ds->GetCellData()->GetScalars()->SetTuple(i, &val);
              }
              else if( atts.GetStatisticType() ==
                  StatisticalTrendsAttributes::Residuals)
              {
                // Slope
                double val =
                  (*(sumXY_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                   (*(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                    *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i))) /
                   (double) (GetNFrames())) /
                  (*(sumX2_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                   (*(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                    *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i))) /
                   (double) (GetNFrames()));
                
                slope_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

                // Average values
                val = *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i)) /
                  (double) GetNFrames();
                sumY_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

                val = *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) /
                  (double) GetNFrames();
                sumX_ds->GetPointData()->GetScalars()->SetTuple(i, &val);

                // Intercept
                val= *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                  *(sumX_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                  *(slope_ds->GetPointData()->GetScalars()->GetTuple(i));

                intercept_ds->GetPointData()->GetScalars()->SetTuple(i, &val);
              }
            }
          }
        }
      }

      // Calculate the sum of squares.
      if( out_ds->GetPointData()->GetScalars() &&
          curr_ds->GetPointData()->GetScalars() )
      {
        unsigned int tPoints = trend_ds[0]->GetNumberOfPoints();
        unsigned int nPoints = curr_ds->GetNumberOfPoints();
        
        if( nPoints && nPoints == tPoints )
        {
          // Traverse all point data
          for( unsigned int i=0; i<nPoints; ++i )
          {
            double val = 0.;

            if( atts.GetStatisticType() ==
                StatisticalTrendsAttributes::Variance ||
                atts.GetStatisticType() ==
                StatisticalTrendsAttributes::StandardDeviation )
            {
              val = *(curr_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i));
            }
            else if( atts.GetStatisticType() ==
                     StatisticalTrendsAttributes::Residuals)
            {
              val = *(curr_ds->GetPointData()->GetScalars()->GetTuple(i)) - 
                (*(slope_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                 timeAxisValue +
                 *(intercept_ds->GetPointData()->GetScalars()->GetTuple(i)));
            }

            double sum = *(out_ds->GetPointData()->GetScalars()->GetTuple(i)) +
              val * val;
            
            out_ds->GetPointData()->GetScalars()->SetTuple(i, &sum);
          }
        }
      }

      if( out_ds->GetCellData()->GetScalars() &&
          curr_ds->GetCellData()->GetScalars() )
      {
        unsigned int tCells = trend_ds[0]->GetNumberOfCells();
        unsigned int nCells = curr_ds->GetNumberOfCells();
        
        if( nCells && nCells == tCells )
        {
          // Traverse all cell data
          for( unsigned int i=0; i<nCells; ++i )
          {
            double val = 0.;

            if( atts.GetStatisticType() ==
                StatisticalTrendsAttributes::Variance ||
                atts.GetStatisticType() ==
                StatisticalTrendsAttributes::StandardDeviation )
            {
              val = *(curr_ds->GetPointData()->GetScalars()->GetTuple(i)) -
                *(sumY_ds->GetPointData()->GetScalars()->GetTuple(i));
            }
            else if( atts.GetStatisticType() ==
                     StatisticalTrendsAttributes::Residuals)
            {
              val = *(curr_ds->GetPointData()->GetScalars()->GetTuple(i)) - 
                (*(slope_ds->GetPointData()->GetScalars()->GetTuple(i)) *
                 timeAxisValue +
                 *(intercept_ds->GetPointData()->GetScalars()->GetTuple(i)));
            }

            double sum = *(out_ds->GetPointData()->GetScalars()->GetTuple(i)) +
              val * val;
            
            out_ds->GetPointData()->GetScalars()->SetTuple(i, &sum);
          }
        }
      }
    }
}



// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::ExecutionSuccessful
//
//  Purpose:
//
//  Programmer: Allen R. Sanderson
//  Creation:   May 07, 2011
//
// ****************************************************************************
bool
avtStatisticalTrendsFilter::ExecutionSuccessful(void)
{
  return (trend_ds.size() != 0);
}


// ****************************************************************************
//  Method: avtStatisticalTrendsFilter::CreateFinalOutput
//
//  Purpose:
//      This method is the mechanism for the base class to tell its derived
//      types that no more time slices are coming and it should put together
//      its final output. This method creates the final output dataset with
//      either all points of the different time slices or the dataset with
//      the particle paths.
//
//  Programmer: Allen R. Sanderson
//  Creation:   May 07, 2011
//
// ****************************************************************************

void
avtStatisticalTrendsFilter::CreateFinalOutput(void)
{

  if( !validTimeAxis )
  {
    std::string msg;

    if( atts.GetTrendAxis() == StatisticalTrendsAttributes::Time )
      msg += "The time ";
    else if( atts.GetTrendAxis() == StatisticalTrendsAttributes::Cycle )
      msg += "The cycle. ";
    else 
      msg += "The time step ";


    msg += std::string("axis values are present but not valid ") +
      std::string("(not in increasing order). ") +
      std::string("The resulting plot may not be correct. ") +
      std::string("Try using another value for the displaying the time axis.");
                       
    avtCallback::IssueWarning(msg.c_str());
  }

  if( trend_ds.size() )
  {
    if( atts.GetStatisticType() == StatisticalTrendsAttributes::Mean ||
        atts.GetStatisticType() == StatisticalTrendsAttributes::Variance ||
        atts.GetStatisticType() == StatisticalTrendsAttributes::StandardDeviation ||
        atts.GetStatisticType() == StatisticalTrendsAttributes::Residuals)
    {
      // Calculate the average or variance.
      if( out_ds->GetPointData()->GetScalars() )
      {
        unsigned int tPoints = out_ds->GetNumberOfPoints();
        
        if( tPoints )
        {
          // Get the next point and update its coordinates if necessary
          vtkDataArray *out_scalars = out_ds->GetPointData()->GetScalars();
            
          // Traverse all point data
          for( unsigned int i=0; i<tPoints; ++i )
          {
            double val = *(out_scalars->GetTuple(i)) / (double) GetNFrames();

            if( atts.GetStatisticType() == StatisticalTrendsAttributes::StandardDeviation )
              val = sqrt( val );

            out_scalars->SetTuple(i, &val);
          }
        }
      }
        
      if( out_ds->GetCellData()->GetScalars() )
      {
        unsigned int tCells = out_ds->GetNumberOfCells();
        
        if( tCells )
        {
          // Get the next point and update its coordinates if necessary
          vtkDataArray *out_scalars = out_ds->GetCellData()->GetScalars();
          
          // Traverse all cell data
          for( unsigned int i=0; i<tCells; ++i )
          {
            double val = *(out_scalars->GetTuple(i)) / (double) GetNFrames();
            out_scalars->SetTuple(i, &val);
          }
        }
      }
    }
    else if( atts.GetStatisticType() == StatisticalTrendsAttributes::Slope)
    {
      vtkDataSet *sumX_ds   = trend_ds[0];
      vtkDataSet *sumY_ds   = trend_ds[1];
      vtkDataSet *sumX2_ds  = trend_ds[2];
      vtkDataSet *sumXY_ds  = trend_ds[3];

      // Calculate the slope via a simple regression model
      // val = (sumXY - (sumX*sumY) / N) / (sumX2 - (sumX*sumX) / N);

      if( sumX_ds->GetPointData()->GetScalars() )
      {
        unsigned int tPoints = sumX_ds->GetNumberOfPoints();
        
        if( tPoints )
        {
          // Get the next point and update its coordinates if necessary
          vtkDataArray *sumX_scalars  = sumX_ds->GetPointData()->GetScalars();
          vtkDataArray *sumY_scalars  = sumY_ds->GetPointData()->GetScalars();
          vtkDataArray *sumX2_scalars = sumX2_ds->GetPointData()->GetScalars();
          vtkDataArray *sumXY_scalars = sumXY_ds->GetPointData()->GetScalars();
          vtkDataArray *out_scalars   = out_ds->GetPointData()->GetScalars();

          // Traverse all point data
          for( unsigned int i=0; i<tPoints; ++i )
          {
            double val = (*(sumXY_scalars->GetTuple(i)) -
                   (*(sumX_scalars->GetTuple(i)) * *(sumY_scalars->GetTuple(i))) / (double) (GetNFrames())) /
              (*(sumX2_scalars->GetTuple(i)) -
               (*(sumX_scalars->GetTuple(i)) * *(sumX_scalars->GetTuple(i))) / (double) (GetNFrames()));

            out_scalars->SetTuple(i, &val);
          }
        }
      }

      if( sumX_ds->GetCellData()->GetScalars() )
      {
        unsigned int tCells = sumX_ds->GetNumberOfCells();
        
        if( tCells )
        {
          // Get the next cell and update its coordinates if necessary
          vtkDataArray *sumX_scalars  = sumX_ds->GetCellData()->GetScalars();
          vtkDataArray *sumY_scalars  = sumY_ds->GetCellData()->GetScalars();
          vtkDataArray *sumX2_scalars = sumX2_ds->GetCellData()->GetScalars();
          vtkDataArray *sumXY_scalars = sumXY_ds->GetCellData()->GetScalars();
          vtkDataArray *out_scalars   = out_ds->GetCellData()->GetScalars();

          // Traverse all cell data
          for( unsigned int i=0; i<tCells; ++i )
          {
            double val = (*(sumXY_scalars->GetTuple(i)) -
                   (*(sumX_scalars->GetTuple(i)) * *(sumY_scalars->GetTuple(i))) / (double) (GetNFrames())) /
              (*(sumX2_scalars->GetTuple(i)) -
               (*(sumX_scalars->GetTuple(i)) * *(sumX_scalars->GetTuple(i))) / (double) (GetNFrames()));

            out_scalars->SetTuple(i, &val);
          }
        }
      }
    }

    std::string newPipelineVariable =
      outVarName + " " + std::string( typeString[(int) atts.GetStatisticType()] );

    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();

    // Set the new data range.
    double range[2] = { FLT_MAX, -FLT_MAX };

    if( out_ds->GetPointData()->GetScalars() )
      out_ds->GetPointData()->GetScalars()->SetName(newPipelineVariable.c_str());
    if( out_ds->GetCellData()->GetScalars() )
      out_ds->GetCellData()->GetScalars()->SetName(newPipelineVariable.c_str());

    GetDataRange(out_ds, range, newPipelineVariable.c_str(), false);

    outAtts.GetThisProcsOriginalDataExtents(newPipelineVariable.c_str())->Set(range);
    outAtts.GetThisProcsActualDataExtents(newPipelineVariable.c_str())->Set(range);

    avtDataTree_p newTree = new avtDataTree(out_ds, 0);
    SetOutputDataTree(newTree);
  }
}
